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Abstract 

To study the dynamics of chemical processes, we often adopt rate equations to observe 
the change in chemical concentrations. However, when the number of the molecules is small, 
the fluctuations cannot be neglected. We often study the effects of fluctuations with the 
help of stochastic differential equations. 

Chemicals are composed of molecules on a microscopic level. In principle, the number of 
molecules must be an integer, which must only change discretely. However, in analysis using 
stochastic differential equations, the fluctuations are regarded as continuous changes. This 
approximation can only be valid if applied to fluctuations that involve a sufhciently large 
number of molecules. In the case of extremely rare chemical species, the actual discreteness 
of the molecules may critically affect the dynamics of the system. 

To elucidate the effects of the discreteness, we study an autocatalytic system consisting 
of several interacting chemical species with a small number of molecules through stochastic 
particle simulations. We found novel states, which were characterized as an extinction of 
molecule species, due to the discrete nature of the molecules. We also observed a strong 
dependence of the chemical concentrations on the size of the system, which was caused by 
transitions to the novel states. 



I. INTRODUCTION: DISCRETE REACTION SYSTEMS 

In nature, there exists a variety of systems that involve chemical reactions. Some are on 
a geographical-scale, while others on a nano-scale. Chemical reactions are an integral part 
of life, including all living forms of life. 

To study the dynamics of reaction systems, we often adopt rate equations in order to 
observe the change in chemical concentrations. In rate equations, we regard the concentra- 
tions as continuous variables; the rate of the reaction as a function of the concentrations. In 
macroscopic systems, there are a vast number of molecules; thus, continuous representations 
are usually applicable. 

When the concentration of a certain chemical is small, fluctuations in the reactions or 
flow can be significant. We often handle such systems with the help of stochastic differential 
equations, in which we regard noise as a continuum description of the fluctuations • Such 
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an approximation is useful when the number of molecules is intermediate. The employment 
of stochastic differential equations led to some important discoveries such as noise-induced 
order noise- induced phase transitions Q, and stochastic resonance 

In stochastic differential equations, still quantities of chemicals are regarded as contin- 
uous variables. Essentially, on a microscopic level, chemicals are composed of molecules. 
The number of molecules should be an integer (0, 1, 2, •■•), which changes discretely. 
Fluctuations are derivatives of discrete stochastic changes; thus, continuum descriptions of 
fluctuations are not always appropriate and can be doubted. For chemicals with a small 
number of molecules of the order of 1, a single molecule is extremely significant; therefore, 
the discreteness in the number is significant. 

Biological cells appear to be a good example. The size of the cells is of the order of 
microns, in which nano-scale "quantum" effects can be ignored. However, in cells, some 
chemicals act at extremely low concentrations of the order of pM or nM. Assuming that the 
typical volume of a cell ranges from 1 to 10'^ /im'^, the concentration of one molecule in the 
cell volume corresponds to 1.7 pM-1.7 nM. It is probable that the molecular numbers of 
some chemicals in a cell are of the order of 1, or sometimes reach 0. 

If such chemicals play only a minor role, we can safely ignore these chemicals. However, 
this is not always the case. In biological systems, chemical species with a small number 
of molecules may critically affect the behavior of the entire system. For example, there 
exist only one or a few copies of genetic molecules such as DNA, which are important 
to characterize the behavior, in each cell. Further, some experiments show that doses of 
particular chemicals at concentrations of the order of pM or fM may alter the behavior of 
the cells (e.g., ^7]). Biological systems also include positive- feedback mechanisms such as 
autocatalytic reactions, which may amplify single molecular changes to a macroscopic level. 
The effects due to small molecular numbers in cells have been noticed only recently, both 
theoretically [113 and experimentally [Tol |. 

At present, we focus on the possible effects of molecular discreteness. To study such 
effects, we should adopt an appropriate method to handle molecular discreteness. Some 
numerical methods to investigate reaction systems that take into account discreteness and 
stochasticity already exist (we briefly review these methods; see Appendix). Among the 
methods, we adopted Gillespie's direct method, which is popular and frequently used. 

Furthermore, some works related to molecular discreteness also exist. For example, Blu- 
menfeld et al. showed that the mass action law may breakdown in a small swtem [ill ll^ . 
Stange et al. studied the synchronization of the turnover cycle of enzymes |l3l Il4l |. 

We regard it important to identify the phenomena for which molecular discreteness is 
essential. 

Through stochastic simulations, we show that discreteness can induce transitions to novel 
states in autocatalytic systems |l5j| . which may affect macroscopic chemical concentrations 

M 
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II. DISCRETENESS-INDUCED TRANSITIONS 
A. Model 

We consider a simple autocatalj^tic network (loop) with k chemicals. We consider Xi 
chemicals and assume 

Xi + — > 2Xi_|_i 

reactions between these chemicals {i = 1, 2, ••• , k; Xk+i = Xi). All the reactions are 
irreversible. 

For the reactor, we assimie a well-stirred container with volume V. The set of Ni, the 
number of X^ molecules, determines the state of the system. The container is in contact 
with a chemical reservoir, in which the concentration of Xi is fixed at Sj. The flow rate of 
Xi between the container and the reservoir is Di, which corresponds to the probability of 
the flowing out of a molecule per time unit^ . 

We can consider the continuum limit as F — > oo. In the continuum limit, the change of 
Xi, the chemical concentration Xi in the container, follows the rate equation 

— fi—lXi—iXi TiXiXi-^\ -\- Di(^Si ^i)? (1) 

where Vi is the rate constant of the reaction Xi + Xi+i —>■ 2Xi+i, and Xq = X^. 

For simplicity, we consider the case with equivalent chemical species, given as — r, 
Di = D, and Sj = s for all i (r, D, s > 0). By this assumption, the rate equation has 
only one attractor: a stable fixed point Xi = s for all i. For any initial condition, each Xi 
converges to s, the fixed point value. Around the fixed point, Xi vibrates with the frequency 
uip = rs/iT. 

If the number of molecules is finite but fairly large, we can estimate the dynamical behavior 
of the system using a Langevin equation, obtained by adding a noise term to the rate 
equation. Each concentration Xi fluctuates and vibrates around the flxed point. An increase 
in the noise (corresponding to a decrease in the number of molecules) merely boosts the 
fluctuation. 

However, when the number of molecules is small, the behavior of the system is completely 
different. First, we investigate the case when k = A, which is the smallest number of species 
to show the novel states described below. 



B. Novel States Induced by the Discreteness 

Subsequently, we investigate the dynamical behavior of the system with a small number 
of molecules. In order to detect the phenomena for which the discreteness of the number of 
molecules is crucial, we employ stochastic simulations. 



^ Di is the diffusion rate across the surface of the container. Here, we choose the flow proportional to V , 
to have a well-defined continuum limit. One might assume the flow proportional to V^/"^, considering the 
area of the surface. By rescaling D, the model can be rewritten into the case with y2/3 for finite 
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Here, we adopt Gillespie's direct method^. The frequency (expected number per unit 
time) 

• of the reaction Xi + X^+i 2X^+1 is Pba = rxiXi+iV = rNiNi+i/V; 

• of the outflow of Xi is Poi = DxiV = DNi] 

• of the inflow of Xi is Pu = DsV . 

In the continuum limit {V — s- oo), these frequencies agree with the rate equation. We 
calculate these frequencies with the current Ni, and stochastically decide when and which 
event will occur next. 

In this case, by an appropriate conversion of V, and t, we can set r and s to be 1 
without loss of the generality (rs/D and sV are the only independent parameters). We 
assume that r = 1 and s = 1 for the purpose of further discussion. The total number of 
molecules in the container, Ntot = X^^ii approximately 4sy on an average. By varying 
we can control the average number of molecules without changing the continuum limit. 

First, we consider the case of a large V , i.e., both the number of molecules in the container 
and flow of molecules between the container and the reservoir are large. As expected, the 
behavior of the system is similar to that of the rate equation with noise. As shown in Fig. 
12 each Ni fluctuates and vibrates around the fixed point value Ni — V (i.e. Xi = 1). This 
is still in the realm of stochastic differential equations. 

However, when V is small, we observe novel states that do not exist in the continuum 
limit. As shown in Fig. |3| continuous vibrations disappear. Furthermore, two chemicals are 
dominant and the other two are mostly extinct {Ni = 0). In Fig. El at t < 520, Ni and A^a 
dominate the system and A2 = A'4 = for the most part. We call such a state the 1-3 rich 
state. Reversely, at t > 520, N2 and A'4 are large and usually A'^i = A^3 = 0. We call this 
state the 2-4 rich state. 

These states appear because of the following reason. In this system, the production of 
Xi molecules requires at least one Xi molecule as a catalyst. If Xi becomes extinct, the 
production of Xi halts. Ni never regains before an Xi molecule flows in. 

In the rate equation (eq. Q), the concentration Xi is a continuous variable, which can be 
an infinitesimal but positive value. The consumption rate of Xi is proportional to Xi itself; 
thus, Xi cannot reach exactly within finite time, even if it can go to asymptotically as 
t 00. 

In fact, the number of molecules must be an integer. Transitions from A'^; = 1 to Ni ~ 
are probabilistic and may happen in finite time. For transitions to occur, it is important 
that the consumption rate of Xi does not converge to at A'i = 1. The average interval of a 
molecule flowing in is 1/DV for each chemical. If D and V are small enough to ensure that 
the inflow interval is longer than the time scale of the reactions, it is likely that the state of 
the system A^^ = 1 drops to A'i = before an Xi molecule enters. 

When Ni reaches 0, A'i+2 is also likely to become 0. For example, if we assume that 
N2 = 0, then A^i is likely to increase because the consumption of Xi halts; N3 is likely 
to decrease because the production of A3 halts. Thus, this results in A^i > N3. When 



^ For systems with a large number of chemical species, the next reaction method would be more suitable. 
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A^i > iVa, the consumption rate of X4 is larger than the production rate of X4; therefore, 
N4 starts to decrease and often reaches 0. When N2 = N4 = 0, all the reactions stop. The 
system stays at N2 = N4 = for a long time as compared with the ordinary time scale of 
the reactions ^/r). This is the 1-3 rich state. 

In the 1-3 rich state, the system alternately switches between A^i > N3 and Ni < N^. We 
consider that the system is in the 1-3 rich state with A^i > N3. One X2 molecule flowing 
in may resume the reactions Xi + X2 ^ 2X2 and X2 + X3 2X3. Generally, the former 
is faster because A^i > N3; hence, N2 is likely to increase. Since N4 = 0, the reactions are 
one-way; iVi decreases and N3 increases. When iVi < iVs, A''2 starts to decrease. Finally, 
when N2 returns to 0, the reactions halt again. The system stays in the 1-3 rich state, until 
Ni reaches 0. In the same manner, the inflow of X4 can switch the system from Ni < N3 
to Ni > in the 1-3 rich state. Consequently, we observe successive switching between 
A^i > N3 and A^i < N^. In the 2-4 rich state, the system switches between N2 > N4 and 
N2 < N4. 

In this manner, even one molecule can switch the system within the 1-3 or 2-4 rich states. 
We name these states "switching states" . 

Now, we investigate some properties of the switching states. We introduce an index 
z = {xi + X3) — {x2 + X4) as a characteristic of the switching states. Around the fixed point 
of the rate equation, z w 0; in the 1-3 rich state, z « 4; in the 2-4 rich state, z w —4. 

The distribution of z is shown in Fig. ^ When V > 128, a single peak appears around 
z — 0, which corresponds to the fixed point. By decreasing V, the peak broadens with 
fluctuations. When V < 32, double peaks appear at z sa ±4, which correspond to the 
switching states. We clearly observe a symmetry-breaking transition between a continuous 
vibration around the fixed point with large V and the switching states with small V. This 
is a discreteness-induced transition (DIT) that occurs with decrease of V, which is not seen 
in continuum descriptions. 

We introduce another index y = {xi + X2) — {x^ -\- X4), which represents the difference 
in concentrations: xi — X3 in the 1-3 rich state and X2 — X4 in the 2-4 rich state. The 
distribution of y is shown in Fig. There are double peaks around y = ±3, which imply 
large imbalances, such as {Ni,Ns) = (3.5y,0.5T^) and (0.5^,3.51^), between A'^i and in 
the 1-3 rich state (as well as the 2-4 rich state). 

C. Single Molecular Switch 

We investigate some properties of the switching states. 

First, we examine how each Ni changes in a switching event. We assume the 1-3 rich 
state with A^i = Nuni, N3 = N^mi, and N2 = N4 = 0. Here, one X2 molecule flows in 
{N2 = I) at t — 0, which starts up the reactions. Assuming that DV is so small that no 
more molecules flow in or out throughout the switching, the total number of molecules A'' 
is conserved at A^ = Nuni + N^ini + 1, and A4 is always 0. We can represent the state with 
two variables, A^i and A^2- 

In this system, only the following types of reactions 

• Xi+X2^ 2X2 and 

• X2 -\- A3 — > 2A3 
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may change Ni; the others never take place. Ni monotonously decreases. When N2 reaches 
0, these reactions halt completely, and the switching is completed. Evidently, Nifin + 
Nsfin = Niini + N^ini + 1, where Nifin and N^fin are the final values of A^i and N^, 
respectively. 

The frequency of the reaction Xi + X2 2X2 is 

fi{Ni,N2) ^ N1N2/V, 

and that of X2 + X3 ^ 2Xs is 

/2(iVi, 7V2) = N2N3/V = N2iN -Ni- N2)/V. 
We obtain the Master equation 

/i(iVi + 1, 7V2 - l)P{N^ + l,N2- l,t) + /2(iVi, iV2 + l)P(A^i, iV2 + 1, i) 
- {/i(A^i, 7V2) + /2(iVi, iV2)} P(iVi, iV2, i) 
^{{Ni + 1){N2 - l)P{Ni + l,N2- l,t) 

+ {N2 + 1)(^ -N1-N2- l)P{Ni,N2 + 1, t) 

~N2iN -N2)P{N,,N2,t)} (2) 

for P{Ni, N2,t), the probability of residence in the state {Ni,N2) at time t. The initial 
condition is P{Ni,m, 1, 0) = 1; otherwise P(iVi, iVa, 0) = 0. 

We can easily follow the Master equation numerically. We investigate the relationship 
between the initial state {NunijNsini) and the final state {Nifim N^fin). 

If the first reaction is X2 + — > 2X3, N2 instantaneously reaches 0, and as a result 
Nsfin — N^ini + 1- The system fails to switch. 

If Niini > N^ini, it is more probable that the first reaction is Xi + X2 2X2- 
Subsequently, the system carries out further reactions. In this case, it is probable that 
Niini ~ Nsfin, i.e., the system swaps A^i and N^, as shown in Fig. Consequently, we 
observe successive switching (as seen in Fig. (JJ . 

When Niini ^ N^ini, the system is likely to reach Nifin — and break the 1-3 rich state. 
A large imbalance between A^i and N3 results in an unstable 1-3 rich state. 

D. Conditions for Switching States 

Now, we investigate the requirements for the transitions to the switching states. The 
rate of residence of the switching states for several D and V is shown in Fig. |S1 For 
approximately DV < 1, we observe the switching states. If DV < 0.1, the system mostly 
stays in the switching states. 

Subsequently, we expect that switching states appear even for large if _D is very small. 
In fact, we observe switching states for V = 10'* as shown in Fig. |^ Furthermore, in this 
case also, a single molecule can induce switching. 

Strictly speaking, if DV is the same, the rate of residence is a little smaller for larger V. 
If V is large, the system takes longer to reach Ni — even if the rate of reaction and the 
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initial concentrations arc the same. Thus, it is less likely to exhibit switching states for the 
same interval of the inflow, 1/DV . 

In general, some reactions are much faster than the inflow. If the number of molecules 
which enter within the time scale of the reactions is of the order of 1 for a certain chemical, 
the reactions may consume all the molecules of the chemical, and the molecular discreteness 
of the chemical becomes significant. In other words, for the effect of the discreteness to 
appear in a system with several processes, it is important that the number of events of a 
process within the time scale of another process is of the order of 1. 

Once the system is in the switching state, it is fairly stable and difficult to escape, espe- 
cially if DV is small. To escape the 1-3 rich state and regain continuous vibration, at least 
one X2 molecule and one molecule should flow in and coexist. It is required that after 
an X2 molecule flows in, an X4 molecule should flow in before N2 returns to 0, or vice versa. 

We put one X2 molecule into the system in the 1-3 rich state at t = 0, and one 
molecule in at t = r. We assume that DV is so small that no more molecules flow in or 
out. Then, we judge whether the system escapes from the 1-3 rich state. In due course, 
the system returns to the 1-3 rich (or sometimes 2-4 rich) state because there is no further 
flow; thus, we should judge at the right moment. Here, if A'i > for all i at i = 8 (i.e., 
waiting about 2.5 times longer than the period of the oscillation around the fixed point), 
we consider that the 1-3 rich state has been interrupted. We measure the probability of 
interruption for various initial conditions and the delay r as shown in Fig. 1101 

The system requires adequate timing of inflow to escape the switching states, which may 
amplify the imbalance between the stability of each state. For example, to escape the 1-3 
rich, TVi > Nj, state, it is required that an X2 molecule flows in, and then an X^ molecule 
flows in with a certain delay r, as shown in Fig. ^| Thus, the frequency of escape from the 
1-3 rich state is approximately proportional to the product of the inflow frequencies of X2 
and X4. If each Di or Si is species-dependent, the stability of the 1-3 and 2-4 rich states 
may strongly depend on DiSiV, the inflow frequency of Xi. 

Furthermore, if at the outset Ni N^, it is difficult for the system to escape from the 
1-3 rich state. In some cases where the parameters are species-dependent, the flows or the 
switching may lead to A^i ~ A'3, which stabilizes the 1-3 rich state. 

These conditions are important to stabilize particular states and affect the macroscopic 
behavior of the system (see Section ^Oj. 

E. Stability and Shape of the Network: fc 7^ 4 case 

To close the section^ we briefly discuss the cases where fc 7^ 4. Figure ITTI shows the time 
series of each Ni for fc = 3, 5, and 6. When k = 6, 1-3-5 rich and 2-4-6 rich states appear. 
However, these states are less stable than the 1-3 and 2-4 rich states for k — A. These states 
collapse when any of the rich chemicals vanish; thus, they are unstable for large k. When k 
is odd (3, 5, • • •), there are no stable states where particular chemicals are extinct. However, 
additional reactions, or a variety of or Si, may stabilize or destabilize the states such that 
it is not always true that loops with odd-fc are unstable. 

The discreteness-induced transitions are not limited in the autocatalytic loop. We apply 
the abovementioned discussions to certain segments of a complicated reaction network with 
a slight modification. 
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III. ALTERATION OF CONCENTRATIONS BY THE DIT 

In the preceding section, we show that the discreteness of the molecules can induce tran- 
sitions to novel "switching" states in autocatalytic systems. 

For the case where fc = 4 with uniform parameters, the 1-3 rich state and the 2-4 rich 
state are equivalent. In due course, the system alternates between the 1-3 rich and 2-4 rich 
states. The long-term averaged concentrations are still the same as the continuum limit 
value, Xi = 1. 

It will be important if macroscopic properties, such as the average concentrations, can be 
altered. We show that the discreteness-induced transitions may alter the long-term averaged 
concentrations. 

A. Model 

Once again, here, we adopt the autocatalytic reaction loop 

+ Xi^i — > 2Xi_|_i 

for the k = 4: species. Now we consider the case where the parameters Di, Sj, or rj are 
species-dependent. In the continuum limit, the concentration Xi is governed by the rate 
equation 

dX i , . / r» \ 

= T{—\X{—\X{ TiXiX-i-^x -\- U{ySi Xij. (^oj 

The rate equation does not contain the volume V; hence, the average concentrations should 

be independent of V. 

As discussed in the preceding section, for the transitions to the switching states to occur, 
it is necessary that the interval of the inflow is longer than the time scale of the reactions. 
In this model, the inflow interval of Xi is ~ 1/DiSiV, and the time scale of the reaction 
Xi + Xi+i 2X^+1 in order to use Xi up is ~ l/viXi+i. If all the chemicals are equivalent, 
the discreteness of all the chemicals equally take effect, and the 1-3 and 2-4 rich states 
coordinately appear at DV « r. 

Now, since the parameters are species-dependent, the effect of discreteness may be dif- 
ferent for each species. For example, assuming that -Di.si < -D2S2, the inflow interval of Xi 
is longer than that of X2. Thus, the discreteness of the inflow of Xi may be significant for 
larger V. 

To demonstrate a possible effect of the discreteness on the macroscopic properties of the 
system, we measure each average concentration afj, sampled over a long enough time to allow 
transitions between the 1-3 and 2-4 rich states, by Gillespie's direct method. Note that every 
Xi does not depend on V in the continuum limit. Generally, in discrete simulations, the effect 
of the discreteness varies with V and alters every Xi. When V is very large, the discreteness 
does not matter and Xi is almost equal to the continuum limit value. In contrast, when V 
is small, the discreteness causes Xi to be very different from the continuum limit. 

We first investigate the case where each Sj is species-dependent (i.e., each inflow rate is 
species-dependent), and each Di and rj are uniform {Di = D, Vi = 1). Later, we briefly 
discuss the case where rj is inhomogeneous. 
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As mentioned in Section III Dl for the effect of the discreteness to appear, it is important 
that the interval of events of a process is of the order of or longer than the time scale 
of another process. As regards the inflow, there are two indices to determine how the 
discreteness appears: 

• The inflow interval, tt-tt, 

• The number of molecules at equilibrium, SiV. 

If the inflow interval, , is longer than the time scale of the reactions, the reactions 
may exhaust the chemical before the chemical enters, and the inflow discreteness becomes 
significant. 

Furthermore, if SiV is smaller than 1, Ni can reach because of the outflow. In such 
cases, the relation between the inflow interval and the outflow time scale is also important. 
The approach time from iV^ = n ^ 1 to A^i = is of the order of -^logn. If the inflow 
interval of the chemical that causes the switching to raise Ni is long enough to allow all Xi 
molecules to flow out, the inflow discreteness may alter the stability of the states drastically. 

From this point of view, we classify the mechanism in cases I, I', and II as follows. 

B. Case I: inflow discreteness and reaction rate 

We start with the simplest case, Si = S3 > S2 = S4. In this case, the rate equation has 
a stable fixed point with \fi : Xi = Si. When V is large, each Xi fluctuates around the fixed 
point, and each average concentration Xi is in accordance with the fixed point value. When 
V is small, Xi depends on V. Fig. E| shows each function of V. The difference 

between xi and X2 increases for small V. 

By decreasing V, first A2 and A'4 reach and the 1-3 rich state appears. To reach 
A^2 = or A'4 = 0, the inflow interval of X2 or X4 should be longer than the time scale 
of the reactions. We set Xi = 0(1); thus, the condition to achieve the 1-3 rich state is 
approximately ^j^, > i; that for the 2-4 rich state is -pJ^, > i. 

If V satisfies jyj^ , 753757 < ^ < lyj^v ' Ds^v ' ^'•^ ^'^'-^ state appears but the 2-4 rich 
state does not. Thus, xi and X3 increase. We actually observed this at V ^ TJT^ ^ shown 
in Fig. ini 

For smaller V that fulfills both the 1-3 and 2-4 rich states, the imbalance between the 1-3 
and 2-4 rich states does not disappear. Once the system is in the 1-3 rich state, adequate 
timing for the X2 and X4 inflow is required to escape the state. Thus, the frequency of 
escape from the 1-3 rich state is approximately proportional to D'^ S2S4V^ , the product of 
the inflow frequencies of X2 and X4. The average residence time in the 1-3 rich state as 
well as the 2-4 rich state is the reciprocal of the escape frequency. The ratio of the average 

residence time in the 1-3 rich state to that in the 2-4 rich state is w . In addition, after 

S2S4 ' 

escaping the switching states, the system tends to reach the 1-3 rich state rather than the 
2-4 rich state because of the biased inflow. Thus, the ratio of the total residence time in the 
1-3 rich state to that of the 2-4 rich state is larger than j^', hence, afi, xz » X2-, X4, even 
for a small difference in Si. 
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C. Case I': imbalance of inflow discreteness 

In Case I, we consider the imbalance between the Xi, X3 pair and the X2, X4 pair. If 
another imbalance exists between the X2 and X4 inflows, the switching induced by these 
chemicals in the 1-3 rich state may be unbalanced. We consider the case where si = S3 > 

S2 > S4. 

In this case, the 1-3 rich state is more stable than the 2-4 rich state, which is identical to 
Case I. In the 1-3 rich state, the system can be switched from iVi > N3 to Ni < N3 by an 
X2 molecule; and from Ni < N3 to Ni > N3 by an X4 molecule. Now, the inflow rate of 
X2 is larger than X4; thus, switching from A^i > to iVi < N3 is more probable than vice 
versa, and the system tends to stay in the iVi < state. Consequently, xi < X3, as shown 
in Fig. El This effect requires switching states. When V is large, xi and X3 are almost the 
same. 



D. Case II: inflow and outflow 

Here, we consider the case where si = S2 > S3 = •S4- In this case also, the rate equation 
has a stable fixed point^. 

When V is small, both the 1-3 and the 2-4 rich states appear. Here, we consider S4V, 
the number of X4 molecules when the concentration of X4 in the container and in the 
reservoir are at equilibrium. If S4V < 1, 7V4 reaches without undergoing any reaction. 
The system takes ~ jj logn time units to reach from N4 = n ^ 1 to N4 = 0. The reaction 
X4 + Xi 2X1 also uses X4 such that iV4 decreases faster if si is large. 

If N4 > 0, the inflow of X3 may switch from N2 > N4 to N2 < N4 and raise A'4 again. 
The inflow interval of X3 is ~ lyj^- the interval is much shorter than the approach time 
to A'4 = 0, the switching maintains A'4 > 0. If the interval is longer, A'4 reaches before 
switching, and the 2-4 rich state is easily destroyed. 

In the 1-3 rich state, the system tends to maintain A'l < N3 because the X2 inflow is 
frequent. However, the Xi inflow is also large enough to maintain A'l > 0. The 1-3 rich 
state retains its stability. 

In conclusion, the 1-3 rich state is more stable than the 2-4 rich state. In the 1-3 rich 
state, A^^i < N3 is preferred, and £3 increases. At this stage, it is possible that £2 ^ £3 
despite the fact that S2 ^ S3, as shown in Fig. H"!] 

E. Amplification by Discreteness 

In summary, the difference in the "extent of discreteness" between chemical species induces 
novel transitions. The "extent of discreteness" depends on V; thus, we observe transitions 
by changing V. The transition reported in Sectionals regarded as a second order transi- 
tion involving symmetry breaking (see Figs. 21and[Sl), while the transition in this section 



^ Generally, if Z) <C rsi, xi, ^ {si + s^)/2 and X2, X4 ^ (s2 + S4)/2. 
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corresponds to the first order transition without symmetry breaking (see Fig. I13|l in terms 
of thermodynamics. 

We classified the mechanism in cases I, I', and II. These mechanisms can be combined. 
For example, we demonstrate the case where si = 0.09, S2 = 3.89, and — S4 ^ 0.01. In 
this case, each Xi shows a three-step change with V, as shown in Fig. 1151 

When V is large, xi, <^ X2, X4, since si + S3 <C S2 + S4. AtV^ 10^, the discreteness of 
the X3 inflow becomes significant, and the 2-4 rich state appears. In the 2-4 rich state, the 
system tends to remain at > -^4 because of the inflow imbalance between Xi and X3, 
as observed in Case I'. Figure [T7I shows the distribution of X2- The major peak corresponds 
to the 2-4 rich, X2 > X4 state for the cases when 128 < V < 512. 

On the other hand, in the 2-4 rich state, the outflow of X4 depresses A'4 toward S4V, 
as observed in Fig. ^| By decreasing V, the imbalance between N2 and A'4 increases 
because the rate of switching, which again raises 7V4, decreases in proportion to V. Finally, 
at F « 10^, the 2-4 rich state loses stability, as seen in Case II. Now, the 1-3 rich state is 
preferred despite the fact that si + S3 ^ S2 + S4. ^3 increases to « 2, which is more than 
30 times as large as that in the continuum limit. 

For extremely small V, the 1-3 rich state is also unstable because A^i and A^3 easily reach 
0. In such a situation, typically only one chemical species is in the container. The system 
is dominated by diffusion, and X2 increases again due to the large S2. 

Note that the chemical that becomes extinct depends not only on the flows but also on 
the reactions. In some cases, we observe smaller x^ for larger S3. 

F. Asymmetric Reaction Case 

It is possible that each varies with the species. In such cases, we can discuss the effect 
of discreteness in a similar way. However, the change of Xi with V is different from the case 
with asymmetric flows. 

For example, we assume that ri = > r2 — ^4 and Vi : s^ = 1. In the continuum limit 
or in the case of large V, £2 = £4 > xi = x^,, as shown in Fig. 1181 In contrast, when V is 
small, Xi K 1. If is very small, such that the total number of molecules is mostly or 1, 
reactions rarely take place. The flow of chemicals dominate the system; thus, Xi ~ Si. 

If both the reactions and the flows are species-dependent, we simply expect the behavior 
to be a combination of the abovementioned cases. Even this simple system can exhibit a 
multi-step change in concentrations along with a change in V. It is not limited to the simple 
reaction loop. In fact, we observe this kind of change in concentrations with a change in 
the system size in randomly connected reaction networks. For a large reaction network with 
multiple time scales of reactions and flows, the discreteness effect may exhibit behavior that 
is more complicated. Our discussion is largely applicable to such cases if we can define the 
time scales appropriately. 

As seen in this paper, the discreteness of molecules can alter the average concentrations. 
When the rates of inflow and/or the reaction are species-dependent, transitions between 
the discreteness-induced states are imbalanced. This may alter the average concentrations 
drastically from those of the continuum limit case. 
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IV. DISCUSSION 

We demonstrated that molecular discreteness may induce transitions to novel states in 
autocatalytic systems, and that may result in an alteration of the macroscopic properties 
such as the average chemical concentrations. 

In biochemical pathways, it is not anomalous that the number of molecules of a chemical 
is of the order of 10^ or less in a cell. There are thousands of protein species, and the total 
number of protein molecules in a cell is not very large. For example, in signal transduction 
pathways, some chemicals work at less than 100 molecules per cell. There exist only one or 
a few copies of genetic molecules such as DNA; furthermore, mRNAs and tRNAs are not 
present in large numbers. Thus, regulation mechanisms involving genes are quite stochastic. 
Molecular discreteness naturally concerns such rare chemicals. 

One of the authors, Kaneko, and Yomo recently provided the "Minority Control conjec- 
ture," which propounds that chemical species with a small number of molecules governs 
the behavior of a replicating system, which is related to the origin of heredity 17, 18, 19] . 
Matsuura et al. experimentally demonstrated that a small number of genetic molecules is 
essential for evolution [23|. Molecular discreteness should be significant for such chemicals, 
and may be relevant to characters of genetic molecules. 

Until now, we have modeled reactions in a well-stirred medium, where only the number 
of molecules is taken into account while determining the behavior. However, if the system is 
not mixed well, we should take into account the diffusion in space. Both the total number of 
molecules and the spatial distributions of the molecules may be significant. From a biological 
point of view, the diffusion in space is also important because the diffusion in cells is not 
always fast as compared with the time scales of the reactions. If the reactions are faster 
than the mixing, we should consider the system as a reaction-diffusion one, with discrete 
molecules diffusing in space. The relation between these time scales will be important, as 
indicated by Mikhailov and Hess 14, 21, 2j|. As regards these time scales, we recent ly fo und 
that the spatial discreteness of molecules within the so-called Kuramoto length 0, l23l |2J| , 
over which a molecule diffuses in its lifetime (lapses before it undergoes reaction), may yield 
novel steady states that are not observed in the reaction-diffusion equations ^5]. There is 
still room for exploration in this field, e.g., pattern formation. 

Our result does not depend on the details of the reaction and may be applicable to systems 
beyond reactions, such as ecosystems or economic systems. The infiow of chemicals in a 
reaction system can be seen as a model of intrusion or evolution in an ecosystem; both 
systems with discrete agents (molecules or individuals), which may become extinct. In 
this regard, our result is relevant to studies of ecosystems, e.g., extinction dynamics with a 
replicator model by Tokita and Yasutomi 26, 27j. The discreteness of agents or operations 
might also be relevant to some economic models, e.g., artificial markets. 

Most mathematical methods that are applied to reaction systems cannot account for the 
discreteness. Although the utility of simulations have become convenient with the progress 
of computer technology, it might be useful if we could construct a theoretical formulation 
applicable to discrete reaction systems. On the other hand, in recent years, major advances 
have been made in the detection of a small number of molecules and fabrication of small 
reactors, which raises our hopes to demonstrate the effect of discreteness experimentally. 

We believe that molecular discreteness is of hidden but real importance with respect to 
biological mechanisms, such as pattern formation, regulation of biochemical pathways, or 
evolution, to be pursued in the future. 
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Appendix: Methods for Simulating Discrete Reaction Systems 

Since the 1970s, several methods have been suggested for simulating discrete reaction 
systems. We briefly review some of these methods: StochSim method, Gillespie's methods, 
and their improved versions. 

These methods are based on chemical master equations. In chemical master equations, 
we define the state of the system as the number of molecules in each chemical; the reaction 
process as a series of transitions between the states^. We consider each event, i.e., a reaction 
event between molecules, inflow or outflow of a molecule, as a transition. They take place 
stochastically with a certain frequency (probability per time unit) determined by the current 
state. 

For the simulations used in this study, we adopted Gillespie's direct method for simplicity'^. 
We also attempted a direct simulation and the next reaction method, and confirmed that 
our result does not depend on the simulation method. 

A. Direct Simulation with Fixed Time Step 

First, let us consider a very simple approximation, which looks similar to the Euler method 
for differential equations. 

In this method, we fix the time step as At. Assuming that the frequency of the event-i is 
Oi, the average number of the event-i for each step is a^At. If ^ aiAt <^ 1, we approximately 
assume that at most one event occurs at each step, and the probability of the event-i is Ai 
(it is possible that no event occurs in the step). We select an event with a random number, 
change the state according to the event, and recalculate each a^. 

B. StochSim Method 

From this simple concept, we can derive some variations. The so-called StochSim method 
[2^ is one such variation. 

The StochSim method adopts random sampling of molecules. For bimolecular reactions, 
we randomly choose two molecules from the system, and decide whether they react or not 



* Gillespie demonstrated tha t the chemical master equations are exact for gas-phase, well-stirred systems 

in thermal-equilibrium |2S|. 
^ The next reaction method would be more suitable if the system contains more chemical species and 

reactions. 
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with certain probability. The method requires three random numbers in total for each 
step. In case there are some single molecular (first-order) reactions, we choose the second 
molecule from the molecules in the system and some pseudo- molecules (dummies). If the 
second molecule is a pseudo-molecule, then we select the single molecular reaction of the 
first molecule, and determine whether it occurs. 

In the StochSim method, At is restricted by the fastest reaction (with the largest reaction 
rate per pair of molecules). If most of the bimolecular pairs do not react with each other, 
or the reaction probability varies with the species, this method may be impractical because 
no reaction occurs in most steps. 

C. Gillespie's Direct Method 

Incidentally, if the system consists of discrete molecules, it is typical that each frequency 
Oi changes only when an event actually occurs. Taking this into account, Gillespie suggested 
two exact simulation methods: the Direct Method 30] and the First Reaction Method [sij . 
In these methods, we do not fix the time step. Instead, we calculate the time lapse until 
the next event. 

In the direct method, first, we consider the total frequency of the events, a = ^ a^. If 
does not change until the next event, the time lapse until the next event, r, is exponentially 
distributed as P{t) = ae"""^ (0 < r). We determine the time lapse t with an exponentially 
distributed random number. 

Subsequently, the probability that the next event is i is a^/a. We determine which event 
occurs with a random number. We then set the time t forward by t, update the state 
according to the event, and recalculate each frequency a^. Iterate the above steps until the 
designated time elapses. 

D. Gillespie's First Reaction Method 

The first reaction method is similar to the direct method. It is based on the fact that 
Ti, the time lapse until the next event-i, is exponentially distributed as P{Ti) = a^e^"''^* 
(0 < Ti). Only the event with minimum Ti actually occurs. We update the state, recalculate 
each ai, and generate all Ti again with the new corresponding a^. 

In the first reaction method, we need as many random numbers each step as the types 
of events. We calculate all r^, choose the earliest, and discard the others. Generally, the 
processor time to generate random numbers is very large; hence, a large amount of time is 
wasted for several types of events. 

E. Next Reaction Method 

To solve this performance problem, Gibson and Bruck proposed a refinement of the first 
reaction method: Next Reaction Method Although, in general, Monte Carlo simula- 
tions require independency of random numbers, they proved a safe way of recycling random 
numbers, which drastically promotes efficiency. 
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In the next reaction method, we store U, the absolute time when the next event-z occurs, 
instead of r^. In the first step, we have to calculate ti = Ti for every i, according to the 
exponential distribution P{Ti) = flje""*'^' (0 < Ti). We choose the event with the smallest 
ti. According to the event, we set the time t forward to U, change the state, and recalculate 
each Oj. 

In steps that follow, we recalculate each U as follows. 

1. As regards the event just executed, we should recalculate ti = t + Ti with the expo- 
nential distribution of Ti, identical to the first step. 

2. As regards other events whose frequency has changed, we should recalculate the 
corresponding ti. For such events, we convert ti as 

^new 

told is the ti before the event and t„etu is that after the event. 

With this conversion, the actual frequency is adjusted from Uoid to a„eun without using 
random numbers^. 

3. As for other events whose frequency has not changed, we do not need to recalculate 
the corresponding ti. 

Subsequently, we choose the smallest t,;, and proceed further. 

If the event executed does not influence Oj, we do not need to recalculate the concerned 
Ui and ti (except for those of the event just executed). Thus, it is useful to manage the 
dependency of on each event. With this intention, we prepare a dependency graph 
that shows which should be updated after event-j (event-i depends on event-j). In a 
large reaction network, such as biochemical pathways, one chemical species can react with 
only a small part of the chemicals in the entire system. In such cases, recalculation is not 
required for irrelevant chemical species; hence, we can accelerate simulations with help of a 
dependency check. 

It is also important to find the smallest ti quickly. For this purpose, we use a heap, a 
binary tree in which each node is larger than or equal to its parent. The root is the smallest 

at any instance'''. 

The next reaction method requires only one random number per event, and executes much 



Note that Gibson and Bruck mathematically proved the vaUdity of the method and did not mention 
numerical errors in their paper. In some cases, repetition of the conversion might be numerically dangerous. 
In case there are various time scales of reactions, incautious coding may result in numerical errors (e.g., 
rare events would not occur forever). 

The heap sort is an 0{NlogN) sorting algorithm. If only one t, has changed in the step, the cost of 
resorting is O(logiV). 
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faster than first reaction method, especially in case of many chemicals and reactions*. 
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FIG. 1: Time series of the concentration Xi at the continuum limit for r = 1, s = 1, and D = 1/32. 
With D > 0, the rate equation has a stable fixed point \/i : Xi = 1. Each Xi converges to the fixed 
point. 




FIG. 2: Time series of the number of molecules Ni for V = 512 and D = 1/256. Each Ni fluctuates 
around the fixed point (as seen in Fig. but does not reach 0. 
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FIG. 3: Time series of Ni for 1/ = 32 and D — 1/256. In this stage, Ni can reach 0, and the switching 
states appear. In the 1-3 rich state, the system successively switches between the A'^i > N3 and 
A''i < states. The interval of switching is much longer than the period of continuous vibration 
(« tt) observed in Figs. 0and|21 Around t = 520, a transition occurs from the 1-3 rich to the 2-4 
rich state. 




FIG. 4: The probability distribution of the index z = {X1+X3) — (12-1-2:4), sampled over 5 x 10® time 
units, {z is actually a discrete value. Here, we show the distribution as a line graph for visibility.) 
D = 1/128. When V is large {V > 256), there appears a single peak around z — 0, corresponding 
to the fixed point state Wi : Xi = 1. For V < 64, the distribution has double peaks around z — ±4. 
The peak z 4 corresponds to the 1-3 rich state, with A'^i -f ^ 4V, N2 ^ N4 ^ 0. z ^ -4 
corresponds to the 2-4 rich state as well. 
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FIG. 5: The probability distribution of the index y = {xi + X2) — {xs + X4), sampled over 5 x 10^ 
time units. D = 1/128. When V is large, there appears a single peak around y = 0, corresponding 
to the fixed point. When V is small and the system is in the switching states, the index y shows an 
imbalance between the rich (non-zero) chemicals. For example, in the 1-3 rich state, y corresponds 
to xi — X3 since N2 — N4 = for the most part. The distribution shows double peaks around 
y — ±3. Assuming that Xi + x^ = i and 12 = 2:4 = in the 1-3 rich state, y = ±3 correspond to 
(xi^xa) — (3.5, 0.5), (0.5, 3.5). Both the chemicals are likely to have a large imbalance. 




FIG. 6: Probability density for the switch from (NijNa) — {Nuni, Nsini) to {Nifi„,N3fi„). We 
assume a switching event triggered by a single X2 molecule in the 1-3 rich state, and we set the initial 
conditions as N2 = 1 and A'4 = 0. Assuming that there is no further flow {D — 0), we numerically 
follow the master equation (eq. (|5J) until t — 50 (sufficiently large to ensure ~ for most cases). 
For each initial condition, we show the transition probabilities of the final states (A'^i/m, A^s/m)- 
V = 32, Nuni + N^ini — 128. The system shows high probabilities around N^fin ~ N^ini + 1 
(immediately terminated) and A^s/m — N^ni (switching; Nuni > N^ini). 
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FIG. 7: Probability for N2 ~ 0, i.e., reactions have already stalled, as a function of t (in other 
words, the cumulative distribution of the time when N2 reaches 0). V = 32, Nuni + Nuni = 128. 
In the case where Nuni < Nzi„i or Nuni ~ A^ami , the probability increases just after the reactions 
start. For such initial conditions, the reactions terminate near the initial state, as shown in Fig. |S| 
If N-s_ini S> Nsini, the probability steeply increases at t ~ 4, which corresponds to the switching. 
The system takes « 4 time units to complete the switching. 




FIG. 8: The rate of residence of the switching states plotted against DV, the inflow frequency, 
sampled over 10*^ {V > 1024), 10^ {32 < V < 1024), and 10* {V < 32) time units. Here, we define 
the 1-3 rich state as continuation of the state in which at least one of N2 or is for 8 time units 
or longer. Thus, the system may contain states with only one or no chemical, especially for small 
V. We observe the switching states for DV < 1. 
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FIG. 9: Time scries of Ni for V = W and D = 10"". For such a small D, we observe the switchiriEi; 
states for relatively large V. In this case also, a single molecule can induce switching. One Xi 
molecule flows in (iVi = 1), propagates to more than 10*, and then returns to 0. 
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FIG. 10: Probability for interruption of the 1-3 rich state as a function of the delay t, sampled 
over 10® times for each condition. We assume the system where N2 = N4 = 0. We inject an X2 
molecule at i = 0, then an X4 molecule at t = r (no further flow). If Vi : M > at t = 8, we 
ascertain that the system has escaped from the 1-3 rich state. V — 32, Nuni + N-^ni ~ ■iV{— 128). 
An initial large imbalance, such as Nuni '■ Nzini = 120 : 8 (15 : 1), makes it easier to escape the 1-3 
rich state. The system is unlikely to escape from the state when Nuni « A^sini or Nuni < Nuni- 
The probability is maximum at t w 2, which approximately corresponds to a half of the period of 
the vibration around the flxed point. 
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FIG. 11: Time series of Ni for fc = 3, 5, 6, with F = 64 and D = 1/256. For the case where fc = 6, 
the transition occurs from the 1-3-5 rich state to the 2-4-6 rich state at t w 10080. For the cases 
where k = 3 and fe = 5, there are no stable states such as the 1-3 rich state when fc = 4 or 1-3-5 
rich state when fc = 6. 
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FIG. 12: The average concentration Xi as a function of V , sampled over 10^ {V > 1024), 10^ 
(32 < V < 1024), and 10* {V < 32) time units (same for Fig. [TIJ. r = 1 and D = 1/128. (a) Case 
I: si = S3 — 1.7, S2 ~ Si — 0.3. When V is small, the 1-3 rich state appears, and the difference 
between xi and X2 increases, (b) Case I': Si — S3 = 1.9, S2 — 0.19, S4 — 0.01. When V is small, 
X2 and Xi decrease, identical to Case I. In this case, the imbalance between xi and X3 appears 
at the same time. (Figures 1121 1141 1151 and 1171 are reproduced from ref. [16|| by permission of the 
pubhsher.) 




FIG. 13: Probability distribution of y = {xi + X2) — (a::3 + a::4), for Case I': si — S3 = 1.9, S2 = 0.19, 
S4 = 0.01, r — 1, and D = 1/128. There is a peak around y — for large V. The difference in Si 
has a slight effect on Xi. For the case where V < 1024, there appears another peak at j/ « —0.8, 
which corresponds to the 1-3 rich, A''i < N3 state. 
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FIG. 14: The average concentration Xi, for Case II: si = S2 = 1.99, ss = S4 = 0.01, r = 1, and 
D = 1/128. For small V, the 1-3 rich state is preferred, and xs increases. 




FIG. 15: The average concentration Xi for si = 0.09, S2 = 3.89, S3 = S4 = 0.01, r = 1, and 
D = 1/64, sampled over 5 x 10'^ {V > 32) or 5 x 10* {V < 32) time units. By decreasing V, first, 
the 2-4 rich state appears, as seen in Case I. Then, the 2-4 rich state becomes unstable and gives 
way to the 1-3 rich state, as seen in Case II. When V is extremely small {V < 0.5), the flow of 
molecules governs the system, and X2 increases again. 
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FIG. 16: Time scries of iV, for si = 0.09, S2 = 3.89, S3 = 54 = 0.01, r = 1, and D = 1/64. (a) 
V = 256. The 2-4 rich state is dominant (Case I'). Inflow of X3 molecules induces switching from 
N2 > N4 to N2 < N4, which prevents N4 from decreasing to 0. (b) V = 32. Now, the X3 inflow 
is rare, which allows A'4 to reach before the switching. Thus, the 2-4 rich state is unstable (Case 
II). 
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FIG. 17: Distribution of X2 for si = 0.09, S2 = 3.89, S3 = S4 = 0.01, r = 1, and i3 = 1/64, sampled 
over 5 x 10^. For large V, a single peak around X2 — 2 appears, which corresponds to the fixed 
point in the contiimum limit. At V ~ 10'*, double peaks appear around X2 ~ 1 and X2 ~ 3, which 
correspond to the 2-4 rich state. By decreasing V, the two peaks spread apart. At F w 10^, the 
skirt of the low-density (left) peak touches X2 = 0, implying that N2 (and N4) is likely to reach 0, 
and thus, the 2-4 rich state loses stability. A peak at X2 = steeply grows by decreasing V further. 
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FIG. 18: The average concentration Xi for Vi : Si = 1 and D = 1/128 with inequivalent reaction 

constants. For small V , the flows of molecules dominate the system. Thus, a^i « 1, which simply 
reflects Sj = 1; this does not depend on how the continuum limit is imbalanced by the reactions, 
(a) ri = rs = 1 and r2 = rA = 0.9. (b) ri = r2 = 2 and rs = r4 = 1. 



